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We introduce a new scenario for heavy ion collisions that could solve the lingering problems 
associated with the so-called HBT puzzle. We postulate that the system starts expansion as the 
perfect quark-gluon fluid but close to freeze-out it splits into clusters, due to a sharp rise of bulk 
viscosity in the vicinity of the hadronization transition. We then argue that the characteristic cluster 
size is determined by the viscosity coefficient and the expansion rate. Typically it is much smaller 
and at most weakly dependent of the total system volume (hence reaction energy and multiplicity). 
These clusters maintain the pre-existing outward-going flow, as a spray of droplets, but develop 
no flow of their own, and hadronize by evaporation. We provide an ansatz for converting the 
hydrodynamic output into clusters. 

PACS numbers: 25.75.-q, 25.75.Dw, 25.75.Nq 
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I. INTRODUCTION 

One of the most unexpected, and as yet unexplained, 
experimental results found at the Relativistic Heavy Ion 
Collider (RHIC) concerns the description of particle in- 
terfcrometry observables [l|. Before RHIC was turned 
on, it was expected that the deconfined matter would be 
a highly viscous, weakly interacting quark gluon plasma 
0. Thus, ideal hydrodynamics would not provide a 
good description of flow observables sensitive to the earl 
stages of the collision, such as azimuthal anisotropy 
as viscous corrections to these observables would be too 
large. Nevertheless, after data was released, hydrody- 
namic simulations offered in fact very successful inter- 
pretation of transverse momentum spectra and their az- 
imuthal anisotropy However, hydrodynamics also 
lead to the prediction, that a clear signature for the phase 
transition would be an increase of the "out" to "side" 
emission radius ratio (referred to as Ro and R^) due to 
longer lifetime of the system, caused by the softening of 
the equation of state in the transition/crossover region 
0. This prediction turned out not to be true Q. 

Measured parameters Ro and Rs are nearly identi- 
cal. Their (positive) difference R^ — R^ is thought 
to correspond — somewhat simplified — to the duration of 
particle emission. Hence, it looks like the fireball emits 
particles almost instantaneously and does not show any 
sign of phase transition or crossover. Hydrodynamics, 
with "reasonable" freeze-out condition (such as a critical 



temperature of 100 MeV or so) can not describe this even 
qualitatively. 

There could be three possible approaches to the HBT 
puzzle. It could be that the system is simply too com- 
plicated, and that once we include all possible improve- 
ments (full 3D calculation, viscosity, hadronic kinetic af- 
terburner, in-medium hadron modifications etc.), every- 
thing will fit. It could also be that we are drastically mis- 
understanding the data, and the HBT puzzle is a symp- 
tom of inapplicability of hydrodynamics to heavy ion col- 
lisions. Finally, it could be that the hydrodynamic ap- 
proach is basically correct, but there is just one clement 
of physics relevant to freeze-out that is fundamentally 
misunderstood. 

The second possibility is unlikely because, in some 
ways, hydrodynamic prescription does fit HBT data. 
Scaling of HBT radii with the multiplicity rapidity den- 
sity (dN/dyY^^ , over a large range of energies 0] is typ- 
ical for an isentropically expanding fluid that suddenly 
breaks apart. The very good description, within param- 
eters compatible with what is needed to describe flow, of 
the azimuthal dependence of HBT radii [1] , also suggests 
that the hydrodynamic framework is a good ansatz for 
describing the matter produced in heavy ion collisions up 
to freeze- out. 

The first possibility appears, however, also problem- 
atic: successful models and/or parametrisations of the 
freeze-out which describe HBT radii are found in the lit- 
erature [1, [13, [O, [3 1 they could provide a way 
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to gain insights into what is missing. However, we feel 
that successful description involves a dynamical descrip- 
tion from initial conditions plus a freeze- out criterion, 
rather than a fit to data with assumptions put in "by 
hand". Such a description is so far lacking^. Further- 
more, the most plausible refinements to hydrodynamics, 
naniely implementation of fully three-dimensional mod- 
els [lal and the addition of a kinetic theory afterburner 
[l^ do not do anything to solve the HBT discrepancy, 
but in fact make it worse, suggesting that the problem 
is not refinements but rather one large missed physical 
effect. 

One such effect discussed so far in the literature within 
a hydrodynamic context is spinodal clustering (T7l.fT8l.fT9l. 
f20j driven by a first order phase transition (some authors 
talk about explosive freeze-out, within the same context 
[2ll,f23,fll,f2l). 

While clustering has been applied to describe HBT 
data [2^, such a description is yet to be made fully ki- 
netic. One reason why such an ansatz has not been ac- 
cepted so far is that lattice strongly indicates that the 
transition at RHIC energies is not of first-order, but 
rather a smooth cross-over, and the critical point ap- 
pears at a considerable chemical potential. Considering 
the rather universal scaling 0, fH, fl^] found in soft ob- 
servables, an alternative which is not dependent on the 
assumption of first order phase transition is desirable. 

In this paper, we propose such a possible alternative 
freeze-out mechanism, where the dynamics explain why 
it can lead to freeze-out considerably different from the 
usually used typical energy density. We argue that it 
can, in fact, be the basis of reconciling hydrodynamics 
and interferometry. Instead of freeze-out happening at a 
critical temperature or energy density, we speculate that 
the system breaks up into fragments, as a result of the 
bulk viscosity sharply rising close to the phase transi- 
tion temperature. This explanation has the virtue that 
it is connected to theoretical features of QGP, namely its 
near-perfect conformal invariance at high (perturbative) 
temperatures, and the existence of a conformal anomaly 
in the non-perturbative regime. 

II. THE BEHAVIOUR OF BULK VISCOSITY 
NEAR Tc 

The bulk viscosity of high temperature strongly in- 
teracting matter has recently been calculated using per- 
turbative QCD [2^, and found to be negligible, both in 
comparison to shear viscosity and w.r.t. its effect on any 
reasonable collective evolution of the system. This is not 
surprising: The QCD Lagrangian, as long as no "heavy" 



Some kinetic models incorporating partonic interactions, such as 
|29ll , manage to reproduce HBT data for certain values of the par- 
ton scattering cross-section. The interpretation of these results 
within a collective picture is however not yet fully understood. 



quarks are present, is nearly conformally invariant [28|. 
Since, within a fluid, the violation of conformal symme- 
try is linearly proportional to a bulk viscosity term [i^ . 
the near conformal invariance of the QCD Lagrangian 
should guarantee that bulk viscosity is nearly zero, in 
the perturbative regime. 

In the hadron gas phase, of course, the numerous scales 
associated with hadrons render conformal invariance a 
bad symmetry, and hence it is natural to expect that 
bulk viscosity is not negligible. 

This is, again, rooted in a fundamental feature of QCD: 
the non-perturbative conformal anomaly, that manifests 
itself in the scale (usually called Aqcd) at which the 
QCD coupling constant stops being small enough for 
the perturbative expansion to make sense. This scale 
coincides with the scale at which confining forces hold 
hadrons together. 

This violation of conformal invariance is not seen 
perturbatively, but should dominate over the pertur- 
batively calculated bulk viscosity as temperature drops 
close enough to the QCD phase transition. 

What happens to bulk viscosity in this regime, where 
hadrons are not yet formed, presumably the matter is 
still deconfined, but conformal symmetry is badly bro- 
ken? While we can not as yet calculate this rigorously, 
there are compelling arguments [13, fH, fl^ that bulk vis- 
cosity rises sharply, or even diverges, close to the phase 
transition temperature. 

Lattice simulations find that (=0 for a conformally 
invariant system), increases rapidly close to T^. Remem- 
bering that the shear (r/) and bulk (C) viscosities roughly 
scale as [H, fs^, fll] 

V ~ Tolasticy (1) 
C ~ Q - ^'s^ Tinclastic?''^ (2) 

where T(ino)ciastic refers to the equilibration timescale of 
(ine)clastic collisions. Assuming Tinoiastic 1/T allows 
to extract the bulk viscosity from the lattice, and yields 
a sharp rise close to Tc. This can be more formally seen 
from finite temperature sum rules in conjunction with 
lattice data [3^. l33l| . 

The rise is, in fact, likely to be considerably sharper 
than suggests. The dependence of Tjneiastic on tem- 
perature can be guessed from the fact that, at Tc, the 
quark condensate {qq) acquires a finite value, and the 
gluon condensate {G^vG^") sharply increases at the 
phase transition. "Kinetically" , therefore, timescales of 
processes that create extra qq and GG pairs should di- 
verge close to the phase transition temperature, by anal- 
ogy with the divergence of the spin correlation length in 
the Ising model close to the phase transition. Numerical 
studies with viscous hydrodynamics and a chiral model 
[sot seem to confirm that the second-order chiral phase 
transition makes the viscosity diverge. 

The sharp rise of bulk viscosity can also be under- 
stood within string kinetics: confinement, microscopi- 
cally, can be thought of as a "string tension" appearing 
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in the potential. Given the small mass of light quarks 
(and hence ease of creating qq pairs), the appearance of 
even small string tension will lead to a small "preferred 
scale" at which strings break. Hence, conformal sym- 
metry should be quickly badly violated right at the dc- 
confinement phase transition. In particular, in a regime 
where the momentum exchange of the average collision 
is more than enough to break the string, the relevant de- 
grees of freedom are still quarks, not mesons, and the 
shear viscosity is still low, a profound change happens: 
each previously elastic collision, that before just diffused 
momentum, becomes inelastic, where the final state has 
less kinetic energy than the initial state. Even if this dif- 
ference (the energy needed to break the string) is low, 
over many collisions, the heat energy would be converted 
into creating more slightly colder, less pressing particles. 
That's exactly the kind of processes that contribute most 
to bulk viscosity [33|. 

These arguments give evidence to the conjecture that, 
close (from above) to Tc, bulk viscosity goes rapidly from 
a negligible value to a value capable of dominating the 
collective evolution of the system. That this transition 
is sharp can be seen by the sharpness of the lattice de- 
confinement transit from lattice flavor correlations stud- 
ies (such as (Ai^AS*) [3^) seem to confirm that, imme- 
diately above Tc, the relevant degrees of freedom be- 
come quasi-particles similar to the asymptotically free 
quarks. It is therefore likely that non-perturbative effects 
(such as the conformal anomaly) go away soon above Tc. 
Conversely, they should appear quite suddenly if T ap- 
proaches Tc from above, in an expanding cooling fluid. 

In the next section, we will show how this picture could 
yield a freeze-out scenario that has the potential to re- 
solve the HBT puzzle. 



III. CLUSTERING AT THE VISCOSITY PEAK 

It has been noted that bulk viscosity could be helpful 
in making Ro and Rs agree with experiment [40| . Of 
course, too large bulk viscosity, in the context of heavy 
ion collisions, would just mean that the approximations 
on which hydrodynamics is based are not accurate (the 
mean free path is not negligible w.r.t. the system's in- 
homogeneities) and an approach not rely ing on assumed 
near-equilibrium might be necessary [20, |30| . 

One guess for the qualitative behaviour of such a sys- 
tem is illustrated in Fig. [T] one just needs to think what 
would happen if a rapidly expanding compressible fluid 
suddenly becomes rather rigid, sticky and resistant to 
further expansion and deformation. ^ If the material 



^ It should be noted that QCD matter is different enough to "ordi- 
nary" matter that concepts such as "solid", "liquid", "sticky" etc. 
are misleading, since these words carry tacit assumptions that 
arc valid in the "everyday" world but potentially badly broken 



spent a considerable amount of time in the low viscos- 
ity phase there will be pre-existing collective flow that 
is pushing the system outwards. The inhomogeneities 
of this flow will rapidly generate strong viscous forces, 
which will tend to decelerate and stop the expansion. 
These forces, by causality, will not be able to quickly 
overcome the pre-existing flow globally^ but more than 
enough to overcome it locally. If the appearance of vis- 
cosity is sharp enough, these forces cannot overcome the 
inertia of the system and it is natural to suppose that 
the system will be rapidly broken apart into small frag- 
ments, each flowing with pre-existing flow, with roughly 
QGP density. 

It should be clearly stated that although this scenario 
shares some superficial similarities with the more usual 
nucleation picture analyzed, e.g., in El, HOj HI] , it is 
physically completely different. 

In [23| , hadronic bubbles form in a steady supercooled 
QGP medium, and the role of viscosity is to dissipate 
the latent heat during bubble growth. The nucleation 
examined here, on the other hand, occurs due to the 
interplay of a suddendly appearing viscosity with the pre- 
existing advective forces. 

The scenario in [20l | requires a first-order phase tran- 
sition, our does not. The nucleation in [23| proceeds via 
creation of critical bubbles whose size is determined only 
by thermodynamical quantities, viscosity influences the 
dynamics of bubble growth but not the critical size, and 
global collective expansion plays no role. Wc require a 
robust pre-existing global expansion, and a sharp peak 
of viscosity at the critical temperature that forces the 
system to disintegrate into fragments which, in principle, 
can depend on both local and global system properties. 

The bubbles in [2^ are made of hadron gas. Ours are 
evaporating droplets of hot QGP. In [2^ , clustering will 
entail an entropy increase, proportional to the latent heat 
and explained by the different entropy density of the two 
phases. In our approach, the formation of clusters should 
quickly kill off d^u^, so entropy generation (9^?/'')^) 



in QCD. In solid state physics, short-range potentials are usu- 
ally dominated by steep "walls" driven by the Pauli exclusion 
principle. Thus, materials with small inter-molecular distances, 
either small- viscosity "good liquids" or large viscosity "solids", 
are almost always highly incompressible, because of the steep 
inter-molecular potential at small average molecular separation. 
This incompressibility is usually assumed in the definition of both 
"liquid" or "solid". In QCD there is no such short-range repul- 
sive potential, the short-range EoS is the conformally invariant 
ideal gas one, and lattice calculations show that compressibil- 
ity of QCD matter is never high even when the bulk viscosity 
rises close to Tc- The highly viscous phase is thus not a "solid" 
or a "glass", since these materials are usually defined as incom- 
pressible, while highly viscous QCD matter can be compressed 
easily if done infinitely slowly. At finite compression/expansion 
speed, on the other hand, the system becomes "solid-like" due 
to the viscosity, an argument central to our subsequent reason- 
ing. Such characteristics are, to our knowledge, highly unusual 
in solid state physics yet well motivated in QCD 
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FIG. 1: (Color online) Fragmentation of the fireball due to sharply increasing bulk viscosity as the temperature decreases. 
Matter which expanded easily before we describe as oil. It suddenly becomes very rigid against expansion (described as honey 
in the figure) and breaks up into fragments. Hadrons evaporate from these fragments. 



during clustering should be negligible. 

Thus, this work and pre-existing clustering (20| should 
not be treated as complementary descriptions of the same 
phenomena, but as competing scenarios to be differenti- 
ated at the theoretical level (arc the conditions for cither 
scenario relevant to heavy ion collisions?) or through 
experimental data (since our clusters are very different 
from those examined in [H [H, HO, HI]) . 

Relativistically, too high a viscosity indicates that the 
system is too far from equilibrium for the Navier-Stokes 
equations to be a good description. At best, higher or- 
der corrections to hydrodynamics, e.g. [3l| become dom- 
inant. At worst (and more likely), the whole expansion 
in flow gradients becomes divergent. The scenario con- 
sidered here represents a guess of how the system could 
evolve after hydrodynamics breaks down as an appropri- 
ate physical description. 

While this guess, with reasonable timescales/cluster 
sizes, should not suffer from the causality violation 
pathologies that affect first order viscous hydrodynamics, 
we do not at the moment see a way to formally assess its 
likelihood from transport theory arguments beyond de- 
riving some quantitative conditions for this scenario to 
be plausible. We shall do so in the rest of this section. 

The first condition is that forces due to bulk viscosity 
must overwhelm advective forces 



ai 



P 



< 1 



(3) 



where p is the pressure and velocity field. If bulk 
viscosity diverges, as argued in this work and (3^ . the 
applicability of this condition is assured. While an ab 



initio lattice extraction involving pure gluons [3J| does 
not have evidence for a divergence, they estimate 



0.5 < -(T = 1.02Tc) < 2 



(4) 



(s is entropy density.) Assuming an inviscid conformal 
Bjorken [14| dynamics until the rise in bulk viscosity, we 
have 



T 

e 



To X 1/3 

r 

T0A4/3 



(5) 
(6) 

(7) 



where e stands for energy density. 

We can then obtain an estimate for ai in terms of the 
initial temperature Tg and thcrmalization timcscale tq 



C 



(8) 



where Tc = t(T = Tc) This might be less than unity 
even in the "worst case scenario" , although, of course, 
the divergence of C would make clustering much more 
plausible. 

The second condition is that the appearance of the 
viscosity divergence is sudden enough for it not to be 
dissipated by hydrodynamic evolution. For a qualitative 
estimate, 



1 dC 
Cd7 



» ^^u^' 



1 dC dr 
CdT VdT 



(9) 
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or the onset of bulk viscosity will just be dissipated 
through hydrodynamic evolution. Once again, a sharp 
divergence of C will ensure that this condition is satis- 
fied. Considering, as a toy model, a Gaussian peak of 
the evolution of ( w.r.t. temperature, and assuming cr^ 
to be the width of the peak, we get 



1 dC _ 2{T - Tc) 1 dr 
CdT ^ 27rcr2 TdT 



(10) 



fitting a cr^ to the output of |32l | and comparing with a 
Bjorken estimate for dT/dr should convince us that this 
criterion is very plausible as T approaches Tc from above. 

The third condition is that terms in the second order 
of the flow gradient are not relevant for the system un- 
der consideration. The effect of these is the emergence 
of a time-scale (the relaxation time, rn), which delays 
the appearence of viscous forces from the built-up of the 
flow gradient. If the divergence in ( is too sharp around 
Tc, or if the relaxation time is too long, the singularity 
in C will have no effect on the dynamics: by the time 
the viscous forces turn on, the system has allready been 
cooled to below Tc and viscosity is not anymore singular. 
If AT is the width of the peak in this condition can 
be quantitatively estimated as 



a2 



Tn dT 



(11) 



Tn is famously difficoult to estimate from quantum field 
theory. The estimate closest to strongly coupled QCD 
we have is provided by calculations in super-symmetric 
Yang-Mills theories^. |35| 



1 -log2 
6nT 



(12) 



Estimating from p32| AT ~ O.lTc, and assuming once 
again a Bjorken-type evolution before the divergence, we 
find 



a2 



10(1 -log 2) 1 

IStT TcTc 



(13) 



since in general, a2 ^ cei the fulfillment of the second 
condition is more likely than the first at all energies. 
The introduction of multi-dimensional expansion (rather 
than the Id case of will for sure increase dfj,u^ and 
dT/dr, lowering ai and raising a2 (Since a2 ^ ai, this 
raises the possibility for the presently considered sce- 
nario) . 

The estimates here, given our very limited understand- 
ing of some key parameters, should only be taken to un- 
derstand that clustering is not outright excluded. Our 



^ Note that this estimate is for a conformally invariant theory, so 
it is in direct contradiction with our scenario, and hence partic- 
ularly unreliable. Still, nothing more realistic is available at the 
moment 



aim in the subsequent sections of the paper is to try 
to quantitatively estimate some phenomenological con- 
sequences of this scenario, and try to connect it to ex- 
perimental data. 



IV. AN ESTIMATE FOR THE CLUSTER SIZE 

The simplest argument is based on the assumption that 
bulk viscosity diverges at the critical point and therefore 
decouples from the problem. The relevant scales are set 
by Aqcd and Tc- Each fragment will have a typical size 
Rc ^QCD (t-'^*' preferred scale of the system) and move 
in the direction determined by its pre-existing collective 
flow field. The typical energy density in fragments is 
about aTc with a = ■k'^vqgp /i^, where vqgp ~ 30 is ef- 
fective number of degrees of freedom in the quark-gluon 
plasma (QGP). For our rough estimates we take the crit- 
ical temperature T^ = 165 MeV and Agcn = 250 MeV. 
Then the typical fragment mass is estimated as 



M (X ^^Rl^vQGpT^c 



2GeV. 



(14) 



Such a cluster (droplet) will decay into about 10 pions 
or a few heavier hadrons. Note that this estimate is 
good for a cluster containing no strangeness or baryon 
number. To handle these, Eq. ([T^ needs to be updated 
to accommodate strangeness and baryon content, per- 
haps using the methods outlined in |41|. Naively, the 
higher energy content of baryon and strangeness rich 
QGP should increase Aiciustor, so that clusters contain- 
ing baryons (strange and non-strange) should also have 
the mass of high-lying baryonic resonances and decay 
into several particles. Note also that T might well be 
considerably larger than the phase transition (let alone 
the chemical freeze-out) temperature. It is simply the 
temperature at which the bulk viscosity starts becom- 
ing strong enough to locally counteract the built-up fiow. 
The large bulk viscosity, collective manifestation of the 
inter-particle confining potential, will prevent these frag- 
ments from expanding further. They should therefore 
be considered Hagedorn-style "fireballs" rather than as 
expanding fluid clumps. Cascading of these fireballs into 
the ground-state hadrons produces the hadrons at "chem- 
ical freeze-out" . 

The fact that the scale suggested here is similar to the 
hadronic scale begs the question of whether this picture 
is significantly different from the "usual" Cooper-Frye 
particle emission picture. The difference is that within 
the Cooper-Frye scenario, the mostly produced particle 
is the "masslcss" pion, while in our scenario only sys- 
tems having "hadronic" ~ ^qcd mass scale are created 
at hadronization. These systems, furthermore, are not 
"particles" (zero temperature states) but rather finite 
temperature fireballs, although it is reasonable for them 
to transform into Hagedorn-type resonances and decay. 

In the presence of collective expansion Aqcd might 
interplay with other scales of the problem set by cxpan- 



6 



sion velocity gradients. Let us use them in an estimate 
of the size of fragments related to the dynamics of the 
expansion. In order to do so it is useful to recall that the 
energy momentum tensor, with vanishing shear viscosity 
but non-vanishing bulk viscosity is 



integral at the critical point, when actual break-up hap- 
pens, then 



Ski, 



2m 



T^' = {e+ p)u^u'' - pg^^'' + C dpUp (g'^" - u'^u'') (15) Taking « r^, wc get finally 



From energy-momentum conservation d^T^^'^ 
then obtain the rate of energy density decrease 

-u^d^,e = ^ — a^it^ . 

e e 



wc 



(16) 



Note that when (^dpU^ ^ p the energy density decreases 
at the same rate as if no work was performed in case 
with vanishing viscosity. For lower rates of the energy 
density decrease the expansion even decelerates. Micro- 
scopically, this is mediated by inter-particle forces which 
hold the system together. It can happen that the inertia 
of the bulk overcomes these forces and the system thus 
fragments. 

In order to obtain a more quantitative estimate of 
droplet size, wc determine it by the balance of deposited 
energy and collective expansion energy. According to the 
definition of viscosity, it determines the amount of energy 
deposited per imit volume and unit time, i.e. 



dv / dTC(a^u^) 



(17) 



where ( is bulk viscosity and collective 4- velocity. For 
simplicity let us assume again the Bjorken [l3 | picture. 
Then — ^/t a-nd the 3- velocity is ~ z/t. If bulk 
viscosity is indeed rapidly divergent at Tc , we can replace 
it with the ^-function 



C(t) = CcTeS (r(T) ~ Te) = CcTe 



dr 
dT 



(18) 

where is a model parameter which should be given 
by deeper theoretical consideration. If we call = 
Tc^It^t^ ^eget 



Edis ~ SL- 



(19) 



where S is the transverse area of the Bjorken cylinder 
and L is the droplet longitudinal size. Wc consider a 
droplet whose center of mass is located at 2; = (though 
this assumption is not really important due to the boost 
invariance of the system) . 

The kinetic energy of droplet 's expansion, which is in 
fact dissipated due to viscosity, is in non-relativistic limit 



Skin = dVeiT)vl 



(20) 



where e(r) is the internal energy density of the fiuid. It 
is of course a function of time but the above expression 
contains only volume integration. Let us evaluate the 



(21) 



(22) 



Notice that in the numerator is actually the inverse 
expansion rate d^u^^. Thus the droplet size squared is 
inversely proportional to the expansion rate. Within this 
scenario the droplet size will grow with the lifetime of 
the hydrodynamic stage (from the initial temperature 
To to Tc), but the growth will generally be slower than 
linear. For our toy model example where the system 
has a conformal equation of state and Boost-invariancc 
{dN/dy - ~ r^), this growth will be - (dN/dyy/^, 
but it is likely to be slower than that when transverse 
expansion is considered. 

Whether the cluster size is indeed only dependent on 
the internal scale of the system Aqc d (Eq. or on an 
interplay between the internal and collective scales (Eq. 
I22p is difficult to determine from first principles, as it 
depends on a quantitative understanding of the details 
of the non-equilibrium evolution around Tc. 

The main point argued in the last section, one that 
docs not depend on these details, is that the sharp rise of 
bulk viscosity could force the system to break up into dis- 
connected fragments, of a scale and lifetime much smaller 
than the size of the system (0(1 GeV)). These clusters 
then flow apart with pre-existing flow velocity and, pre- 
sumably, decay by Hagedorn cascading. In the next three 
sections we shall examine the effect this kind of freeze-out 
has on heavy ion phenomenology. 



V. PHENOMENOLOGY OF CLUSTERING 

While the HBT puzzle is our main experimental moti- 
vation for introducing a qualitatively new freeze-out sce- 
nario, several observables, aside from particle interfcrom- 
etry, could imply clustering. In this section, we point out 
a list of such phenomena. In each of these, the evidence 
for clustering is by no means overwhelming, and alter- 
native explanations for each of these phenomena exist. 
Nevertheless, it is worthwhile to point these phenomena 
out individually as candidates for contact between the 
model presented here and experimental data. 

• As pointed out in [46] , a highly viscous but hydro- 
dynamic evolution is constrained experimentally 
by multiplicity measurements. The dependance of 
multiplicity on centrality has been shown to be 
well described through exclusively initial conditions 
(Glauber model, or, at high energies, the Color 
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Glass Condensate). Since expansion of a viscous 
fluid generates entropy at the rate [l^ 



(23) 



too much viscosity at any stage during the hy- 
drodynamic evolution would spoil the agreement 
between experimentally observed multiplicity and 
ansatze based on initial conditions. This is a po- 
tential problem of all attempts of solving the HBT 
puzzle through viscous but hydrodynamic evolu- 
tion [43,[43 

Viscosity-driven clustering would not have this 
problem, since within the cluster all relative mo- 
tion is very quickly killed. Thus, while ^ might 
diverge, df^u^ would vanish. Since entropy produc- 
tion rate is quadratically proportional to the latter, 
we would expect the entropy content of the system 
to not be significantly changed during the cluster- 
ing and freeze-out phase. 

• The very fact that a "single freeze-out model" 
mi works much better than naively expected 
in describing soft observables in heavy ion collisions 
suggest that "something" is decreasing hadronic in- 
teractions after chemical freeze-out below the ex- 
pected rate. Clustering of the system into smaller 
sub-systems that decay after a finite time, dur- 
ing which they flow out with the pre-existing flow, 
would have just such an effect. 

• The ower-abundance of certain resonances (S*, A, 
S*) with respect to even chemical freeze-out ex- 
pectations [i^, [4^ could also be nicely explained 
in terms of clustering. Clusters can be considered 
as highly excited Hagedorn tower resonances. It 
is therefore natural to suppose that they could de- 
cay through a cascade down the Hagedorn "tower" , 
and hence through the production of resonances. 
Hence, ratios such as S*/A and S*/2 would be 
correspondingly enhanced. Stable particles should 
still be well described by the statistical model: 
the final hadron abundance will be a collection of 
a large number of microcanonically decaying fire- 
balls, carrying grand-canonically distributed en- 
ergy and quantum numbers. 

• The scaling of px fluctuations provides direct ev- 
idence that particles are emitted from clusters, 
containing a small ('^ 5) number of particles in- 
dependently of collision energy or ccntrality \Ed^ . 
The undcr-prediction, by the equilibrium statisti- 
cal model, of fluctuations of ratios such as K/tt [51| 
compounds this evidence, since cluster emission 
would enhance fluctuations of multiplicity yields 
and ratios. The forward-backward multiplicity cor- 
relations [5^ and angular correlations in Cu-|-Cu 
collisions at RHIC [53 also indicate the presence 
of clusters. 



• Clustering into fragments of size about 1 GeV 
could provide an explanation for the invariant mass 
systematics of the inverse slopes Toff observed at 
SPS and RHIC energies, both for stable particles 
(Fig. 1 of 01) and electromagnetic resonance de- 
cays (p p'^l^^ , Fig. 1 left panel of [i^). All in- 
verse slopes for particles less massive than roughly 
1 GeV seem to rise with mass, as expected, approx- 
imately, for "blue-shifted" thermal emission where 
the inverse slope combines temperature and flow 
(Toff ~ T+{vt)M). For masses larger than 1 GeV, 
however, Toff is nearly independent of mass. Clus- 
tering, if all clusters have a mass of about 1 GeV, 
could provide a natural explanation for this obser- 
vation: Below the invariant mass of 1 GeV, parti- 
cles are predominantly emitted from a single clus- 
ter, and hence maintain the memory of that clus- 
ter's flow. Above a mass of 1 GeV, however, parti- 
cles have to be emitted either before cluster forma- 
tion, or through cluster fusion, and hence there is 
occasion for the flow to be "forgotten" . 



VI. LOOKING FOR CLUSTERS IN HBT 

We start by noting [1, [sj] that, in the out-side- long 
coordinate system, HBT radii are directly related to the 
system's spacetime correlation tensor^ 



Rl{K) = ((A.x.f) 



(24) 



+ (^) (,A„») (25, 



and, for pairs of particles having zero net longitudinal 
momentum 



(26) 



where the k vector is the sum of the two momenta (the 

first element, ko, is ~ V m? 4- fc^). Averaging is done 
using the emission function 

/ A{x)S{x,K)d'^x 



{A) {!<) 



J S{x, K)d^x 



(27) 



Ri is straight-forwardly related to the longitudinal length 
of the fireball. A correct treatment of deviations from 
boost-invariance should therefore also contribute to an 
improvement of current discrepancy with experimental 
data^. 



Here I ("long") is the z direction (parallel to the beam), o ("out") 
is the direction of the pair momentum, and s ("side") is the cross 
product of the previous two. 
^ Please refer to 01, on the current status of hydro-experiment 
HBT comparisons 
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As remarked in the Ro ~ Rs result is not easy 
to reconcile with naive hydrodynamics plus a straight- 
forward (critical temperature) emission because: 

• The higher the energy, the longer the emission time, 
the larger is the expected discrepancy between Ro 
and Rs- If the system starts close to the mixed 
phase, the timcscalc of freezing out should be longer 
still due to the softest point in the equation of state. 
Hence, a generic prediction from Eqs. ([24|) and ([25|l 
is that Ro/Rs > 1, broadly increases with energy, 
and exhibits a peak when the energy density is such 
that the system starts within, or slightly above the 
mixed phase. This is in direct contrast with exper- 
imental data, where Ro/Rs ~ 1 is a feature at all 
reaction energies. 

• Generally in a hydrodynamic model the (AxAt) 
correlation is negative, since particles on the outer 
side are the first to freeze-out. This increases 
Ro/Rs further (cf. eq (|25p ). Time dilation due to 
transverse flow does not help enough, as calcula- 
tions show. 

It is immediately apparent that clustering can help solv- 
ing both of these problems. 

• Cluster size, density and decay timescale, is ap- 
proximately independent of cither reaction energy 
or centrality, as can be deduced from Eq. ((22|l . 
Hence, the near energy independence of the (com- 
paratively short) emission timescale, and hence of 
Ro/Rs, should be recovered. 

• If the decay products do not interact (or do not in- 
teract much) after cluster decay, it can also be seen 
that (AxAt) can indeed be positive: outward clus- 
ters are moving faster, resulting in time dilation. 
This effect can be offset by time dilation of clus- 
ter decay by increasing the temperature at which 
clusters form, or by increasing cluster size. 

Recovering the linear scaling of the radii with 
{dN/dyY^^{^ -/Vciustcrs) [3] J while maintaining the cor- 
rect Ro/Rs is also possible if the clusters decay when 
their distance w.r.t. each other is still comparable to their 
intrinsic size. 

Quantitative calculations are necessary before deter- 
mining whether these constraints can be satisfied. The 
technical details of how to perform such calculations, 
from a hydrodynamic code output with a critical tem- 
perature and cluster size, are outlined in the Appendix. 
Hydrodynamics output is needed to specify the cluster 
fiow array and emission array I]*^, (defined in Eq. 

The bulk-viscosity-driven freeze-out adds another pa- 
rameter to ab initio HBT calculations: in addition to crit- 
ical temperature/energy density, we now have the cluster 
size. To see whether this helps solving the HBT problem, 
output from hydrodynamics with a high (T ~ Tc) freeze- 
out temperature should be fragmented into clusters with 



a certain distribution in size, which then produce hadrons 
according to the prescription in the Appendix. 

If this ansatz, and a reasonable mean/variance do re- 
produce the observed Ro,Rs from a realistic hydrody- 
namics output, it would provide a strong motivation for 
looking for clusters in cvent-by-cvcnt physics. Cluster- 
driven symmetry breaking should also lead to distinctive 
signatures in the multipole expansion of the correlation 
function [sH] 

VII. DISCUSSION AND CONCLUSIONS 

We have described a mechanism to generate fragments 
that is solidly grounded in QCD, and does not require 
a first order phase transition. Hence, it is possible 
that hadronization is governed by this mechanism in all 
regimes where an approximately locally thermalized de- 
confined system is produced. 

Potentially, this mechanism can solve the HBT prob- 
lem by adding a further "free parameter" to the system: 
the cluster size. Using the methods described in Section 
IVIl it is possible to see whether a given cluster distri- 
bution, matched to the hydrodynamic output with the 
freeze-out criterion tuned to cluster formation, could re- 
produce the measured HBT radii. 

Future work in this direction includes both a quantita- 
tive comparison between HBT data and the model (with 
a proper hydro input), as well as signatures for clustering 
in event-by-event physics. 
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APPENDIX A: EMISSION FUNCTION FOR 
CLUSTERS 

The HBT emission function of a fiuid breaking up into 
identical clusters which then decay should be given by a 
sum of cluster emission functions 

5(x,p)=^5,(x-4,p) (Al) 

i 

At the cluster rest frame we suppose, in accordance with 
the ansatz described in Section IIIIl that Si is a simple 
Gaussian with no further structure or fiow. This cluster 
decays, also via a Gaussian function after a time r ~ 
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A 



QCD 



after formation, a scale also similar to it's radius. 



Normalizing to the number of particles per cluster, and 
in the Boltzmann approximation 



1 1 



(27r)3 



g-£;7Tg-(i'"W+y'=+z'")/(2r^) 



(A2) 

Collective (cluster) velocity and the hypersurface on 
which clusters are generated are expressed as 



f coshyLi coshj/Ti 
sinh yTi cos 9i 
sinh yTi sin di 
\ sinhj/LiCOshyTi 



tfi cosh yLi 
r-fi cos 0i 
Tfi sin 9i 
tfi sinhj/Li 



(A3) 

Note that E here is used in a somewhat different way than 
in the context of hydrodynamics. In hydrodynamics, 



is defined as the space-time locus of particle emission. In 
the clustering scenario, it describes the space-time locus 
of cluster formation. Note, in this respect, that and 
are not fields, but rather arrays of four- vectors, incor- 
porating a finite set of cluster flow velocities and emission 
coordinates. 

Putting everything together, in the lab frame 



where 



1 l^__E'/Tg-(^'=-x°)(x;3-xo,3)A^A^/(2r^) 



(27r)3 T 



and the Lorentz matrix is 



(A4) 



(A5) 



A'' = 



/ 7TC0shyL 
7t wt cosh cos 

7twt cosh j/l sin 6 

y 7t sinh 



7tWt cosh j/L cos 7tWt coshj/i sin6' 

-J COS^ r} COS 6 sin 6 r. 



cos sin 6 ^ 



cosh tanh 
i;^+tanh^ yL 



1 



. Vj. sm P r) 
ii^+tanh^ yz, 
cos ^ tanh yx. /3 
u^+tanh^ yz, ^ 



7t sinh \ 

tjt cos tanh y^ /□ 
ii^+tanh^ yz, 
sin ^ tanh y^ 



ii^-f tanh^ yz, ''^ 



1 



tanh i/i. 



and 



and of course 



/3 = 7t cosh j/L - 1 



i?' = coshyL7T(£' — Ptvt cos (A0) — pL tanh^i) 

where is the relative angle between the direction of the emitted particle and the motion of the cluster. 
We note that the emission function of each cluster is in the Gaussian form 



Si{x,p) ^ cxp 



where 



(A6) 



(A7) 



(A8) 



(A9) 



We also need a "map" of clusters, giving us the flow and freeze-out time of cluster i. Assuming boost invariance and 
"global" azimuthal symmetry, as well as small cluster size w.r.t. system size, and even distribution of clusters, we get 



^S^{x,p)=^Yl ^d'^t ^'^nkl^^Oi 



ikl 



k=l 1=1 m=~Ny 



XIX! X Si{tf{Tk,y„i),Xf{rk,0i),yf{rk,9i),Zf{Tk,y,n),l3T{rk)) (AlO) 



k=l 1=1 



-Ny 



r 



where 



rk 



'TV, 



2Trl 



2my„ 



(All) 

(A12) 
(A13) 



and the parameters Prirk) and tf{rk) need to be ob- 
tained by "freezing out" a hydrodynamic simulation with 
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the appropriate temperature (the temperature, in the 
QCD phase, where bulk viscosity becomes dominant). 

We can estimate the number of clusters in each direc- 
tion by requiring each "side" of the cluster to be approx- 
imately of length T (note that Nq depends on k) 



(A14) 
(A15) 
(A16) 



Nr 






T 




2TTrk 




T 


Ny 


sin 


= Tfc 



We note that the number of clusters times the number of 
particles per cluster is equal to the total multiplicity so 



{N)^ 



rfnax^max sinh(2y)T 3 



Nr. 



4n 



T^m^TK2 { ^ 



(27r)3 

as expected from the statistical model 



(A17) 
(A18) 



To parametrize tf,(3T.L = tanh(j/j'^i) we use the usual 



Vf 



Tk cosh y,„ 
Tk cos(6';) 
rk sm{0i) 
Tk sinh Um 



(A19) 
(A20) 
(A21) 
(A22) 



(3t and Tk can be obtained through a hydrodynamic cal- 
culation, assuming freeze-out occurs when T Tc- 

Throughout, we shall use the usual "out-side-long" co- 
ordinate system, where the x axis points "outwards" , the 
y axis "sideways" and the z axis "longitudinally" 

In the mass-shell projection (fc^g^ = 0) 




(A23) 



w^/cp = cosh j/i cosh y^fco — sinh j/T cos 0fco — sinh ?/i cosh j/T^L (A24) 

M^f/u = I cosh j/L cosh yTT- — sinh yr cos 9]qo — sinh yx sin Oqs 
\ ko J 

cosh y i cosh yx- sinh y l cosh yr ] ( A25) 

ko J 



and, finally 



Up,XQ = T + tf coshyx ~ r sinh yx (A26) 
q,,xi; = q^j:'' + Tq^u^' (A27) 

ko , ki 
ko ko 



' ko ki 

g^S^ = ( -^qo + —qL ) tf cosh(2/L) - rqo coa{6) - rq^ SYa{6) - qLtf sinh(?/i) . (A28) 



1. A Gaussian approximation estimate 

We use cylindrical symmetry of the emitting function and define 6* = to be in the direction of k (the average 
particle pair momentum). We immediately remember the standard Gaussian integration formula 

*,(fc) = / d'xS,{x^k) = e-<'^v/^^^ (A29) 

It is not difficult to prove that 

/ x^'Six, q)d^x = f Si{x, q)d^x = ^m'^iik) (A30) 

I I 

I d^xx^x^S{x,q) = Yl (-2^'^&^ +<i^o;*'(fc)) (A31) 
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Putting everything together we have 

El 



(A33) 



which can be used in conjunction with Eqs. ([SS]) . ([24|) . ((26|) to calculate Ro,s,i- Note that because of clustering the 
emission function can not he either cylindrically symmetric or boost-invariant, and acquires ojf-diagonal terms (sTIS^]: 



Rl = {Ax,Ax2) ~ ^ {AxoAx2) 

Kq 



(AxoAxa) 



Rl = {AxoAxs) 



(Ax.Ax,) - ^ (AxoAxi) ■ 

fco 



ko 



(AxqAxo) 



while in the cylindrically symmetric case i?os = Rsi = Roi = 0. 
Explicitly, the determinant of B^u is given by [53 | 



and 



cosh(^/T)'*(-10 - 2cosh(2;/L) + cosh(2 (y^. - yr)) + 2cosh(2yT) + cosh(2 (t/l + Vt))? 

64 



3h(yT)^ (cosh(yL)^ cosh(yT)^ + sinh(j/L)' 



-Boo 

i?oi = 2 cos(0) cosh(yL)^ cosh(yT)^ sinh(yT 
Bii 



^fl^2 \2 - / x2 . . ^2 , cos(6')^cosh(yi)%inh(2/i)^sinh(yT)^ 

cos(y) cosh(2/L) cosh(yT) smh(yT) H — 



cos(6')^ cosh(yi)'* sin(6')^ sinh(yT)^ 



(1 + cosh(yL) cosh(yT)) 

cos(6')^ ( — l + cosh(yL) cosh(?;T)) sinh(j/T)^ 
sinh(yT)^ + tanh(yL)^ 



B22 — 



2__^,_,_ n2 _,__,^^2 \2 , cosh(2/L)^ sin(6')^ sinh(?/L)^ sinh(?/T)^ 



(1 + cosh(?;L) cosh(t/T))' 



B33 

B12 
B13 
B23 



(1 + cosh(yL) cosh(yT)) 
cosh(yi)^ cosh(yT)^ sin(6')^ sinh(yT) 

^cos(6')^ cosh(yL)'' sin{9)'^ sinh(yT)"'^ 

(1 + cosh(yL) cosh(yT))^ 
2 + 6 cosh(2 + cosh(2 (j/l - yr)) - 2 cosh(2 yr) + cosh(2 {yL + yr)) 



2 cosh(yi)^ cosh(y7-)^ sin(0) sinh(yx) 
cosh(yT)^ sinh(2j/L) 

cosh(?/i)^ (3 + cosh(2 sin(2 0) sinh(yT)^ 



+ 



-l + cosh(?/L) cosh(?/T)) sin(0)^ sinh(?/T)^ 
sinh(j/T)^ + tanh(yL)^ 



cos(6') (3- 


4 

h cosh(2 2/T)) sinh(2yL) 


sinh(yT) 


sin (61) (3 ^ 


4 

- cosh(2yT)) sinh(2?/i) 


sinh(?;T) 


4 



(A34) 
(A35) 

(A36) 



(A37) 

(A38) 
(A39) 

(A40) 



(A41) 

(A42) 

(A43) 
(A44) 

(A45) 
(A46) 
(A47) 



2. Calculation of the full correlation function 

The full correlation function is given by [s^ 

C{k,q) = 1 + 





Y.tStik,q) 


2 







(A48) 
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We can again use the standard formulae regarding Fourier transforms of Gaussians, where 

\BIJ 



(A49) 



where 



G'(<7") = exp 



(A50) 



Note that, as a result of interference between clusters, the correlation coefficient will pick up an interference pattern 
(which will disappear if the spatial size r varies considerably from cluster to cluster). The FT of the whole system is, 
up to a factor of 2t{2'k)'^ (that cancels). 



E lul I ^'(g^) + E ' /._. .._ ^ COS {q,{x^,^-xU)Gl{qnG,n{qn (A51) 



the bottom of the correlation is, up to the same factor. 



E*^ l^+a*^ 






-K+^JfcMcoshf^"^"""''"^ 



2T 



I IB™ 



(A52) 



If clusters all have the same size, the cosine term would 
give a characteristic oscillating pattern. However, for a 
general cluster distribution such terms should in general 
interfere destructively. Still, the form of the tail should be 
considerably different from the Gaussian approximation, 
especially in the high go,; tail. 

The explicit form of the exponent in the Fourier trans- 
form, q>'{B'-^^y^q'' , which goes into Eq. (|A50P is of 



course a closed formula, but too long to be included here. 

Generalizing this formalism to non-central, non boost- 
invariant collisions is straight-forward. It is also possible 
to generalize this approach to a distribution of cluster 
sizes, by updating Eq. (|A2[) with r ^ and Eq. (|A1[) 
with J2 ^ fi'''i)^ij where /(r,) is the cluster prob- 
ability distribution. For realistic cluster distributions, 
however, Monte Carlo methods might prove necessary. 
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